! SK INLET/OUTLET 1D CHARACTERISTIC NONREFLECTING BCS  

subroutine usdef1( un, u,  ub, ucb, ucsb  &
                     , pn,p, pb,pcb, pcsb &
                     , ron,rho,rhob,rhocb,rhocsb &
                     ,ds1,esm,dt,dtmin,epsilon2,p00,tout,ubg1,finlet1)
	!IMPLICIT REAL *8(a-h,o-z)
use Data_cfoam, only: epsMaximumPrinciple, nonReflectOutletBC
implicit none

    real *8 esm(3)
	real *8  un(3), u(3),  ub(3), uf(3), ucb(3), ucsb(3), ucf(3), ucsf(3)
	real *8 pn,p, pb, pf, pcb, pcsb, pcf, pcsf
	real *8 rhon,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf
	real *8 ds1,ds2,dt,sound,rofon,dtmin,tout
	real *8 ut(3,9)
! SK EXTRA PARAMETERS SPECIAL FOR THE OUTLET
    real *8 uf1,uf2,uf3,pout
    integer :: finlet1

! SK EXTRA VARIABLES
    real *8 qb,qcb,qcsb,rf,rcf,rcsf,rwellb
	real *8 qwellb,qmaxb,qminb,rwellf,rmaxf,rminf
	real *8 gl,gr,rnr,rnl,qnl,qnr,char1f,char1b
	real *8 smachl,smachr,smach,gn,ubg1(3),ubg(3)
! SK

! END EXTRA PARMS
    real *8 ei,z1,vtmp,dttmp
    real *8 char3,q,r,qn,rn,u1x,ron,vf,wf,vb,wb
    real *8 char3f,soundf,gf,sf,qf,qcf,qcsf,char2f,qwell,qmaxf,qminf
    real *8 char3cf
    real *8 char3csf
    real *8 char3b,soundb,gb,sb,rb,rcb,rcsb,char2b,rwell,rmaxb,rminb
    real *8 char3cb
    real *8 char3csb
    real *8 charc1,charc2,charc3
    real *8 a(3,3)
    real *8 epsilon
    real *8 es(3)
    real *8 defb
    real *8 alfp,dtmint,gam,gam2,gamx,gamy,p00,roout,s,scb,scsb,sn,cv
    integer j
    
    external defb
    real *8 epsilon2

    real(8) :: sn1,sn2, ut12,ut13,ut22,ut23, pn1,pn2, ron1,ron2
    real(8) :: alfa

    epsilon=epsilon2
    !if (pcsb.lt.0.95d5) epsilon=0.2
    !if (pcsb.gt.1.15d5) epsilon=0.2
    gam=1.4d0
    gamx=2.0d0*dsqrt(gam)/(gam-1.0d0)
    gamy=((gam-1.0d0)/(2.0d0*gam))
    gam2=2.0d0*gam
    es=esm/dsqrt(esm(1)**2+esm(2)**2+esm(3)**2)
!	p00=101325.0d0 
	cv=717.65d0


	call coor(es(1),es(2),es(3),a)
! 1-u 2-ub 3-uf 4-ucb 5-ucsb 6-ucf 7-ucsf
    call  transform(u,ut(1,1),a)
    call  transform(ub,ut(1,2),a)


    call  transform(ucb,ut(1,4),a)
    call  transform(ucsb,ut(1,5),a)

!    ubg(1)=ubg1
!    ubg(2)=0.0
!    ubg(3)=0.0

    ubg(1:3)=ubg1(1:3)

! SK TRANSFORM THE BACKGROUND VELOCITY
    call  transform(ubg,ut(1,3),a)
!if(finlet1.eq.1.0d0)then
!    ut(1:3,3) = -ubg(1:3)
!end if
! END OF TRANSFORM BCG VELOCITY 

    char3 =ut(1,1)
    char3b=ut(1,2)
 
    char3cb=ut(1,4) 
    char3csb=ut(1,5)


! SK SWITCHES BETWEEN OUTLET/INLET -- SAME AS IN THE PREVIOUS VERSION////////////////////////
!	alfp=0.0
!!	alfp=1.0
    if(nonReflectOutletBC)then
        alfp=0d0
        if(char3csb+ut(1,3).lt.0)alfp=1.0d0
    !    call GetAlfa(char3csb+ut(1,3),alfp)
	    if(finlet1.eq.1) alfp=1.0d0
	else
	    alfp=1d0	    
	end if

    pout=p00*alfp+pcsb*(1.0d0-alfp) ! nonreflecting boundary
    roout=p00/(tout*cv*(Gam-1.0d0))

!    if(ut(1,3).gt.0)alfp=0.0d0

! SK //////////////////////////////////////////////

!    char3 =ut(1,1)*(1.0d0-alfp) + ut(1,3)*alfp
!    char3b = ut(1,2)*(1.0d0-alfp) + ut(1,3)*alfp
!    char3cb = ut(1,4)*(1.0d0-alfp) + ut(1,3)*alfp
!    char3csb = ut(1,5)*(1.0d0-alfp) + ut(1,3)*alfp


    s=dlog(p)-gam*dlog(rho)
    sb=dlog(pb)-gam*dlog(rhob)

    scb=dlog(pcb)-gam*dlog(rhocb)
    scsb=dlog(pcsb)-gam*dlog(rhocsb)

    gb=gamx*dexp(scsb/gam2)
    soundb=dsqrt(gam*pcsb/rhocsb)

! SK COMPUTE THE REFERENCE FORWARD STATE

   uf1=ut(1,5)*(1.0d0-alfp)+ut(1,3)*alfp
   uf2=ut(2,5)*(1.0d0-alfp)+ut(2,3)*alfp
   uf3=ut(3,5)*(1.0d0-alfp)+ut(3,3)*alfp
   rhof=rhocsb*(1.0d0-alfp)+roout*alfp
   pf=pcsb*(1.0d0-alfp)+pout*alfp

!   uf1=ut(1,3) - (1.0d0-alfp) * 2d0 * ut(1,3)
!   uf2=ut(2,3) - (1.0d0-alfp) * 2d0 * ut(2,3)
!   uf3=ut(3,3) - (1.0d0-alfp) * 2d0 * ut(3,3)
!   rhof= roout
!   pf= pout

   sf=dlog(pf)-gam*dlog(rhof)
   gf=gamx*dexp(sf/gam2)
   soundf=dsqrt(gam*pf/rhof)
   char3f=uf1


! END FORWARD STATE


! SK THEN JUST COPY THE RELEVANT PART FROM THE MAIN USDEF ROUTINE (WITH ALL EPS=THE SAME)

    qf=char3f-gf*pf**gamy
    qnl=qf

!SK  + 4 TERMS IN THE OTHER DIRECTION  (B)
    q=char3-gb*p**gamy      
    qb=char3b-gb*pb**gamy
    qcb=char3cb-gb*pcb**gamy
    qcsb=char3csb-gb*pcsb**gamy
! SK + EXTRAPOLATION AND CORRECTION IN THE OTHER DIRECTION (B)

    char2b=char3csb-soundb	
	qwellb=(qcsb-qcb)/(0.5d0*dt)+char2b*(q-qb)/ds1
	qmaxb=dmax1(q,qcb,qb)+dt*qwellb + epsMaximumPrinciple
	qminb=dmin1(q,qcb,qb)+dt*qwellb - epsMaximumPrinciple
	qnr=(2d0*qcsb-qb*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(qnr>qmaxb) qnr=qmaxb
	if(qnr<qminb) qnr=qminb

! SK

	r=char3+gb*p**gamy
	rb=char3b+gb*pb**gamy
	rcb=char3cb+gb*pcb**gamy
	rcsb=char3csb+gb*pcsb**gamy

	char1b=char3csb+soundb
	rwellb=(rcsb-rcb)/(0.5d0*dt)+char1b*(r-rb)/ds1	
	rmaxb=dmax1(r,rcb,rb)+dt*rwellb + epsMaximumPrinciple
	rminb=dmin1(r,rcb,rb)+dt*rwellb - epsMaximumPrinciple

    rnr=(2d0*rcsb-rb*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(rnr>rmaxb) rnr=rmaxb
	if(rnr<rminb) rnr=rminb
!SK +4 TERMS IN THE OTHER DIRECTION (F)
        rf=char3f+gf*pf**gamy
        rnl=rf


!  SK INVARIANTS DEFINITION END


! SK ALL CHAR'S AT THE RHS BELOW MUST BE REFERRED TO THE CENTRES 
!    sound=soundf+soundb        ! not used VSemiletov
!    charc3=char3f+char3csb
!    charc1=char3+sound   
!    charc2=char3-sound         ! not use, VSemiletov



! SK DEFINITION //////////////////////////////////////
    gr=gb
    gl=gf
!    char3f=char3csf  ! not use, VSemiletov
!    char3b=char3csb  ! not use, VSemiletov
! END DEFINITION//////////////////////////////////////

! SK THE BASE CASE: SUBSONIC LEFT/RIGHT GOING WAVES//////////////////////

    gf=gl
    gb=gr
    rn=rnr
    qn=qnl ! NB: VS: slip place
    
!    rn = rnr
!!    !gr = gb
!    qn = qnl
!!    !gq = gf
!    u1x=(qn+rn) * 0.5d0

    u1x=qn+gf*(rn-qn)/(gf+gb) ! NB: VS commented
!    u1x = u1x*(1.0d0-alfp) + ut(1,3)*alfp
    call GetAlfa(u1x, alfa)
    if(alfa>=1.0d0)then
        sn=sf
        ut(1,8)=0.0d0
        ut(2,8)=uf2
        ut(3,8)=uf3
        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
        ron=(pn/exp(sn))**(1.0d0/gam)
    else if(alfa<=0.0d0)then
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
        ut(1,8)=0.0d0
        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1d0/gam)
    else
        sn1=sf
        ut12=uf2
        ut13=uf3
        pn1=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
        ron1=(pn1/exp(sn1))**(1.0d0/gam)
        
        sn2=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
        ut22=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
        ut23=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
        pn2=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
        ron2=(pn2/exp(sn2))**(1d0/gam)
        
        sn = alfa * sn1 + (1.0d0-alfa) * sn2
        ut(1,8)=0.0d0
        ut(2,8) = alfa * ut12 + (1.0d0-alfa) * ut22
        ut(3,8) = alfa * ut13 + (1.0d0-alfa) * ut23
        pn = alfa * pn1 + (1.0d0-alfa) * pn2
        ron = alfa * ron1 + (1.0d0-alfa) * ron2
    end if
    
!    if(u1x<=0.0) then   ! if(charc3<=0.0) then ! SK changed, 2011-05-17
!        sn=sf
!        ut(1,8)=0.0d0
!        ut(2,8)=uf2
!        ut(3,8)=uf3
!        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
!        ron=(pn/exp(sn))**(1.0d0/gam)
!        !if (k.eq.5051) write(14,*) 'f',pn,ron
!        !if (k.eq.5051) write(14,*) 'fqr',q,r
!    else 
!        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
!        ut(1,8)=0.0d0
!        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
!        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
!        pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
!        ron=(pn/exp(sn))**(1d0/gam)
!        !if (k.eq.5051) write(14,*) 'b',pn,ron
!        !if (k.eq.5051) write(14,*) 'bqr',q,r
!    end if

! SK END OF SUBSONIC ////////////////////////////////////////////

! SK RIGHT-GOING SUPERSONIC/////////////////////////////////////////////


    if(char1b>0.0 .and. char2b>0.0) then
        gb=gr
        gf=gb
        rn=rnr
        qn=qnr

        u1x=0.5d0*(rn+qn)
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
        ut(1,8)=0.0
        ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
        ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
        pn=((rn-qn)/2.0d0/gb)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1d0/gam)
    end if


! SK END OF R/G SUPERSONIC ////////////////////////////////

! SK LEFT-GOING SUPERSONIC/////////////////////////////////////////////

    if(char1b<=0.0 .and. char2b<0.0) then
        gb=gl
        gf=gb
        rn=rnl
        qn=qnl

        u1x=0.5d0*(rn+qn)
        sn=sf
        ut(1,8)=0.0
        ut(2,8)=uf2
        ut(3,8)=uf3
        pn=((rn-qn)/2.0d0/gf)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1.0d0/gam)
    end if

! SK END OF L/G SUPERSONIC ////////////////////////////////
   

    call transformr(ut(1:3,8),ut(1:3,9),a)
	do j=1,3    
        un(j)=ut(j,9)+u1x*es(j)
	end do

	vtmp=dsqrt(un(1)**2+un(2)**2+un(3)**2)
    sound=dsqrt(gam*pn/ron)
    dtmint=ds1/(sound+vtmp)
	if (dtmint.lt.dtmin) dtmin=dtmint ! definition of Time step

    return      
end subroutine usdef1